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inentation of networked, real time, three-dimensional battlefield simulations on low 
cost, commercially available graphics workstations. The most recent system, NPSNET, 
has improved in functionality to such an extent, that it is considered a low cost version 
of the Defense Advanced Research Project Agency's (DARPA) SIMNET system. In 
order to reach that level, it was necessary to economize in certain areas of the code so 
that real time performance occurred at an acceptable level. One of those areas was in 
aircraft dynamics. However, with “off-the-shelf” computers becoming faster and 
cheaper, real-time and realistic dynamics are no longer an expensive option. The real- 
istic behavior can now be enhanced through the incorporation of an aerodynamic mod- 
el. To accomplish this task, a prototype flight simulator was buut that is capable of 
simulating numerous types of aircraft simultaneously within a virtual world. Beside be- 
ing easily incorporated into NPSNET, such a simulator will also provide the base func- 
tionality for the creation of a general purpose aerodynamic simulator that is particularly 
useful to aerodynamic students for graphically analyzing differing aircraft's stability 
and control characteristics. This system is designed for use on a Silicon Graphics work- 


station and uses the GL libraries. 


il 


II. 


IIL. 


4 fDe 
C74 (le 
TABLE OF CONTENTS 

INTRODUGTION .. 22 2eceecccccccccciesden acd | 
A. MOTIVATION winnie eee... 1 
B. FOCUS vicccccisscassscseenseoesenis eee 2 
1. Complexity of the Aerodynamic Model-..........:..-7e-e) aus 2 

2. An Appropriate Orientation Model?.........22.0-.-...3eeee 3 

3. Defining Multiple Aircraft Aerodynamics Characteristics... S 

C. SUMMARY OF CHAPVERS <iipc ir eee ses een 4 
AERODYNAMIC MODEL. i recio piers ccs > 
A. INTRODUCTION ...0.0....205 isda eee 5 
B. COORDINATE SYSTEM o...0.06....05. Sees... cus tect 6 
C. “DEFINITION OF TERMS (teense 7 
D. MODEL DESCRIPTION ....sssssssccssssssssecssssssssesssssssntensssssseceensssseseevensseese 1 
J. Aerodynamic Forces and Moments ..............:cc:nnsess--.s..-+.:.2k 2 

2. Calculating Linear and Angular Acceleration .....,.............-..........ee 14 

3. Calculating Linear and Angular Velocity 223)5..... sun 15 

4. Updating Position isiiiic.c0 orem oiler 16 
METHODS OF ORIENTATION fie. ee: .: ieee iene I 
A. INTRODUCTION 2.0, cic ek science ne <cosee ete Ly 
B. SURVEY OF MEVHODS 2 18 
1, Buler Method... ccc verncesevsncseesne ess ccsa ee ete rare ee en 18 

2. Direction COSines .ceiicccQieiesns.s<alsenceeee ee teem Oe <a 20) 

3. Quaternion Method. core recsee eee ee 2) 

a. Defining Quaternions :....-. 2.0. cisco. eee 22 

b. Quatemmion. Multiplication... 9:2 eee 2s 

c. Quaternions in Tenns of Direction Gosines................ ae ee 24 

d. Derivative of the Unit Quatemion = 29.7... eee a 

e. Euler Angles from Quatemions |... 22. ee 26 


Oe ee Noto INE) DISADVANTAGES 00.0... c.ccciccsciesesccctesccsevesceessnescoes 26 
rete we OREN ATION MODEL cic. ......cccceceeescecccsscsnseeesecseeeseceseers 28 

See ire ole Le WR IA SL) AOE epee ccercsececsce he sceccccaeceascascessvcececnsscacscncececscesnenes 30 
Pm, eso eer es mere NY Oneida es cae is decade adlan sila ceteseosessvevconcnsnaees: 30 

fe ee Ee eye te) AO ACTIN UE PIE oo coeec icc cccccsssesscccdevdevevsasess coceossesdeesooess 2 

eed he erty Fea) TA Be IE eee vs stevca es acovsaepcreaveneun es taascccdsvensecsceeens 31 
Ameren) ce Welty TF de EIN PU) tis erence seedety nese casesec sis sbdocvasesidesecoesaceseosssescseees a5 

aL ol een es PE hh occ. assccsdsanes-d\.erseddadasnaydeeaasrtceseeeeeconds OO 

Poe CT OMe ey ALO WM AMON WUGUET COLITO! ...y.51.cc.5000c0ssn0s venworen cr vsccavscvenwoansess 34 

Pao eR OM@TINON-PICL@PED AIRCRAFT i. o.cccccccecceccccestBetesesseeeesnee SD 
Fens E) OiTeme al OT) YIN CO IB ts chide sess cc eeecseccssacesssacnnessuescteens 36 

V. SOMO SIONS ANID PUMORE WORK Siac. ...c tte ons ss ecicesdeseeseeceeenes a7 
PEN DIX Ay AERODYNAMIC MODEL IMPLEMENTATION .................0.0sccccereeess 38 
er PENDIX B: QUATERNION AND MATRIX FUNCTIONS ....000.......cccccccccscsessessenees 4] 
eee ae VIPLEMENTATION INTO NPSNET oocociccccc.....cceeenccoutess ses cccessdneenuaans 44 
APPENDIX D: PHOTOGRAPHS OF IMPLEMENTATION i SI 
MEO EME eer IN Cb reenter tare rte te te sscrcc cnt Mettcc re ccccccaccscsccccccescccncneséscccssseccccudescccceness 54 
Se etme eure evsibtecitey AL) IN EUS ere sce ck ccuceh ereeet eet stte ace eee ree e eevee se cnansessssesedesssceeeas DD 


Rroupeeta: 


Figure 2.1 
Piguye 222 
Rigure 25 
Figure 2.4 
Figuie 2 
Figure 2.6 
Figure 2.7 
Figure 2.8 


Figure 3.1 
Figure 3.2 
Figure 3.3 
Figure 3.4 
Figure 3.5 


Figure 4.1 
Figure 4.2 
Figure 4.3 
Figure 4.4 


LIST OF FIGURES 


Flow of Data and Summary of@hesis*@hapters .....) ce... 4 
Basic Aerodynamic Model ....22zeees sss: 6 
World and Body Coordinate Systems ..-.............c9ee 7 
Terms Defined within the World Coordinate System ...............c.ccccceeeeee ees 8 
Terms Defined within the Aircraft Body Coordinate System ................... 8 
Notation with Respect to Body AxeS 32... 22ers 2 
Terminology Defining Aircraft Controls .....0.05................- 9 
Aircraft Dimensional Specifications:.................-s9s0e.)..---..-- 10 
Aircraft Stability Coefficients .3,2.--.-<....:::-:)-ser ee 2 
Euler Attitude Angle Rotation ...ccssscessssecsessoessceesscertsereesseeereaseese 19 
Direction Cosines in Terms of Euler Angles ........2..0......---.2- 20 
Quatemion Orientation. .5.22......c-c0ssasss sos d eee 22 
Four Parameter Method (2295. <..c.eesscccs tee eae sn 24 
Efficiency Comparison of Euler and Quaternion Methods ..................... oe 
Prototype Flight Simulator Basic Structure 272................... ee 30 
Data Input Records 2. ic51....220tsto..00ssntte .sshonseenneseneste=seees <=» = ee 32 
Aircraft Flight Data Structure .......:..5.:0...2eueceeece-eseese-ceee-2o eer 39 
Algorithm for Computing Time Step.......22.92--- 36 


vi 


I. INTRODUCTION 


A. MOTIVATION 


The current state of the art in sunulation technology has provided today ‘s military with 
many extremely valuable training experiences that could not have been obtained elsewhere 
and, as a result, has greatly increased survivability and readiness. From flight simulators, 
which allow a pilot to explore the edge of the flight envelope without endangering crew or 
multi-million dollar assets, to battlefield sumulators, which allow entire fighting divisions 
to practice command and control without having to incur the enormous costs of running a 
full blown field exercise, computer simulation has become a vay of doing business within 


the military. 


One simulation system designed by the Defense Advanced Research Projects Agency 
(DARPA) is the Simulation Networking (SIMNET) [Thorpe,1987]. SIMNET 1s a 
networked battlefield simulator that allows multiple use interaction on the battlefield at 
many different levels. Vehicle simulators, such as tanks and aircraft, connect to the network 
and become part of a three dimensional world. On the lowest level, this system allows the 
vehicle drivers a chance to study defensive and offensive tactics. On higher levels, 
Company, Battalion and Regiment Commanders can experiment with the movement of 


forces, thereby learning what will and won’t work at their particular level. 


At the Naval Postgraduate School (NPS), an effort to develop a SIMNET type system 
based on commercially available, general purpose, graphic workstations, has been active 
for a number of years. This system, NPSNET, consists of Silicon Graphics workstations 
attached to a local area Etheret [Zyda.et al.,1992]. Eventually, NPSNET will become a 
node on the SIMNET network. 


Over the past few years, as workstations became available that were both faster and 
cheaper. capabilities, that were once too computationally expensive to incorporate, have 
become feasible. One such feature 1s realistic vehicle dynamics. As a result it is now 
possible to enhance the behavior of the system by the incorporation of realistic vehicle 
dynamics. A ground vehicle dynamics model currently in development at NPS will provide 
a greatly improved system for ground vehicle interaction with the enviornment. Because 
aircraft dynamics differs in many respects to ground vehicle dynamics, a Separate dynamic 
and orientation model for air vehicles 1s required. This work provides such a model for 


incorporation into NPSNET. 


Most of the research and development of aerodynamic simulators has involved either 
the development of large multi-million single aircraft platforms utilizing specialized 
computer hardware to increase processing speed, or the development highly proprietary 
flight simulation programs for use on a personal computers. As a result, a gap occurred 
where a general purpose flight simulation program operating on low cost graphics 
workstations, is needed. A system that allows the flight characteristics of any aircraft be 
modeled on a low cost, general purpose, graphical workstation, would be of great use to 


aerodynamic students in studying the stability and control of different aircraft. 
B. FOCUS 


There are several issues to be addressed when incorporating an aerodynamic model 
into a computer simulation. The complexity of the aerodynamic model, which orientation 
model to use and how aircraft data should be represented in the system are the critical issues 


and are center of focus in this work. 
J. Complexity of the Aerodynamic Model 


Of primary importance is that the complexity of the model fit the objective of the 


simulation. A complete aerodynamic model that includes fully articulated control surfaces 


i) 


and atuflow divergence patterns over the aircraft would seriously affect real time 
performance on any computer. Such models are usually computed in non-real time on super 
computers and are not appropriate for use on low cost graphics workstations. On the other 
hand, modeling the dynamics of an aircraft kinematically so that the aircraft's velocity and 
orientation are a linear result of control input, does not reproduce the nuances of aircraft 
motion and response that a user of a flight sunulation would expect. The aerodynamic 
model’s complexity must provide as much realism as possible without reducing the frame 
rate below an acceptable level. Therefore, the focus in creating the aerodynamic model ts 
to create an aerodynamic model of sufficient complexity that provides as much 


aerodynaniic realism as possible without adversely effecting real-time performance. 
2. An Appropriate Orientation Model 


The choice of rotational model is fundamentally limited to either an Euler angle 
or quaternion approach. This has been the subject of heated debate among computer 
scientists as to which rotation model is the most appropriate [Goldiez and Lin,199] |. 
Basically, either model can be used to direct orientation, but, depending on the type of 
vehicle being modeled, one method has certain advantages over the other. The primary 
focus in defining an appropriate orientation model is in determining which provides the 


right tool for the job [Shoemake,1985]. 
3. Defining Multiple Aircraft Aerodynamics Characteristics 


A general purpose flight simulator, such as is to be incorporated into NPSNET, 
should be capable of simulating an unlimited number of different aircraft types. Because 
each type of aircraft exhibits its own specific aerodynamics and handling characteristics, it 
is desirable to change these characteristics depending on which type of aircraft is to be 
piloted. One solution is to base the aerodynamic model on the aircraft’s stability 
coefficients, inertial coefficients, and airframe specifications, all of which are generally 


available by consulting aerodynamic stability and control textbooks. Stability coefficients 


provide a very accurate model of aircraft flight behavior. However, care must be taken 
when modeling some of the newer generation fighters. To umprove maneuverability, these 
aircraft have been designed aerodynamically unstable. Their stability coefficients reflect 


this instability. 
C. SUMMARY OF CHAPTERS 


Chapter II covers those considerations for incorporating an aerodynamic model into 
NPSNET. A detailed aerodynamic model utilizing stability coefficients 1s presented. 
Chapter III investigates the difference between Euler angle and quatemion representations 
for defining rotations. The advantages and disadvantages of both methods are discussed 
and a model is presented for integration into NPSNET. Chapter JV provides 
unplementation details of the aerodynamic and orientation model into the prototype flight 
simulator. Chapter V covers conclusions and lists requirements and suggestions for future 


work in this area. 


Environmental Aerodynaniic Forces, 
Factors Model Moments 
(Chapter 2) (Chapter 2) (Chapter 2) 


Orientation 
Model 
(Chapter 3) 


NPSNET Implementation 
(Chapter 4) 





Figure 1.1: Flow of Data and Summary of Thesis Chapters 


H. AERODYNAMIC MODEL 


A. INTRODUCTION 


In modeling the dynamic behavior of an aircraft within the framework of a computer 
based flight sumulation, one usually begins by creating a mathematical model. This model 
represents the relationship between an aircraft and its interaction with the air mass in which 
it operates. Depending on the complexity of the model, additional forces, such as those 
associated with propulsion and environment, and complex systems, such as flight controls 


and weapons, are included [Rolfe,1986]. 


How detailed this model becomes ts dependent on how it will be utilized. Generally, 
aerodynamic engineers will break the model down into components and study the results 
of detailed, but partial simulations. In complex flight simulators, such as those utilized to 
train pilots and aircrew, the model becomes more complex due to the increased number of 
systems that must be simulated and the response speed that 1s required for proper feedback 
to the pilot. As a result, the model loses fidelity, with the amount of detail limited by 


hardware and real-time constraints. 


Flight dynamics in NPSNET requires a mathematical model similar to that found in 
complex flight simulators. The framework for such a model is presented in Figure 2.1. By 
adjusting the detail of those functions which surround the basic mathematical model, a 
simulation 1s built which can be made more or less detailed based on the speed of the 


workstation on which it is implemented. 


FORCES Input 


Aerodynamic Flight 
MATH MODEL Simulation 
Environmental 
Propulsive 
Airframe 
Specs 


Figure 2.1: Basic Aerodynamic Model 





B. COORDINATE SYSTEM 


Coordinate systems and the method in which they are described vary to a great extent 
on the application and the preference of the user. In aircraft sunulations, coordinate systems 
fall into two broad classes, “body” coordinates and “earth” or “inertial” coordinates 
[Rolfe.1986]. Body coordinates have their origin based at the center of gravity and 
continually move with the aircraft. Inertial coordinates, on the other hand, are defined with 
respect to the earth and have their origin positioned at some suitable location such as the 
center of the simulated world. Other coordinate systems exist, based on parameters such as 
the flight path and angle of attack. However, the aerodynamic model presented in this paper 
is based on geometric body and world coordinate systems. In general, all aerodynamic 
forces, accelerations and velocities are calculated in the body coordinate system first, and 
then converted to the world coordinate system prior to updating an aircraft’s position and 


attitude. 


Figure 2.2 shows the generally accepted convention for labeling of the axes in the two 
coordinate systems found in most aerodynamic textbooks [Anderson,1989]. Body 


coordinates are defined with the origin at the center of gravity (CG), the x axis along the 


fuselage pointing out the nose of the aircraft, the y axis along the wing-line pointing out the 
right wing, and the z axis pointing out the bottom of the plane. World coordinates are 
defined with the origin based at a fixed point on the ground, the x axis pointing North, the 
y axis pointing East, and the z axis pointing down. Because of its limited effect, the 


curvature of the earth is ignored. 


Body Coordinates World Coordinates 


North 


Down 





C. DEFINITION OF TERMS 


In a dynamics model, velocities, accelerations and forces are described in both world 
and body coordinate systems. Without an explicit description of the variables used to 


describe the model, confusion can arise. Most terms described in this paper refer to the 


geometric body axes. However if a reference is made to the world coordinate system the 


subscript “w” is used. See Figure 2.3 for terms defined in world coordinates. 


Aircraft location in world coordinates (feet) 


Aircraft velocity in world coordinates (ft/sec) 





Figure 2.3: Terms Defined within the World Coordinate System 


Terms without the “w” subscript relate to body axes and include linear and angular 
velocities, accelerations, forces, moments, angle of attack and sideslip angle. Figure 2.4 
provides a description of these terms. Note that the direction of angular accelerations and 
velocities and moment terms are defined using the night hand rule around their respective 


axis (Figure 2.5). 
Linear velocity along X, Y, and Z body axes (ft/sec) 
Angular velocity along X,Y, and Z body axes (rad/sec) 
Resultant velocity Vector Creve 


Wind velocity across tail of aircraft 


Linear acceleration (ft/sec”) 


Angular acceleration (rad/sec’) 


Forces acting on aircraft 
Moments about the X, Y, and Z axes 


Angle of attack[tan! (W/U)] 


Sideslip [tan'! (V/U)] 





Figure 2.4: Terms Defined within the Aircraft Body Coordinate System 


The aircraft control surfaces such as elevator, ailerons, and rudder are defined as a 
rotation in radians around their respective hinge points on the aircraft. When a control 
surface is flush with the aircraft, the angle of deflection 1s zero. See Figure 2.6 for a further 


description. 





Figure 2.5: Notation with Respect to Body Axes 


elevator deflection positive down (radians) 
a positive de produces a positive lift and a negative pitch moment 


auleron deflection positive left (radians) 
a positive da produces a negative roll moment 


positive nose left a positive or produces (radians) 
a positive sideforce and a negative yaw moment 





Figure 2.6: Terminology Defining Aircraft Controls 


Within the aerodynamic model, the particular aircraft being modeled is characterized 
by certain dimensional characteristics. A description of these terms is included in Figure 


Le 

surface area of wing (ft’) 
wing span (ft) 

cord length (ft) 

weight (ibs) 

roll inertia (slug-ft’) 


= oO TD 


~ 
~ 


pitch inertia (slug-ft*) 


yaw inertia (slug-ft?) 


N 
N 


I 
In 
I 
I 


roll-yaw cross inertia (slug-ft*) 


~ 
N 





Figure 2.7: Aircraft Dimensional Specifications 


The model presented in this paper utilizes English Standard Units of measurement. 


Therefore, the following units are utilized: 


Distance: . feet 
Force: pounds 
Mass: slugs 
Angles: radians 
Temperature: Kelvin 
Time: seconds 


In NPSNET distances are measured by meters. A conversion can be made into the 
metric system by either changing all the units of measurement to standard metric units prior 
to calculating any result, or converting the final result to meters prior to updating position 
using a simple conversion factor of 0.3048 meters/feet. In this paper the latter method is 


utilized. 


D. MODEL DESCRIPTION 


As indicated in Figure 2.1, the mathematical model presented takes forces, control 
inputs and aircraft specifications as inputs, generating linear and angular velocities in 
aircraft body coordinates as outputs. Based on a classical representation of linear 
aerodynamics and utilizing the total force Equations (2.18 - 2.20), all the forces associated 
with lift and drag are calculated utilizing aerodynamic stability derivatives [Roskam,1979]. 
Stability derivatives, first used by Bryan over a half-century ago, assume that all 
aerodynamic forces and moments can be expressed as a function of the instantaneous value 
of the perturbation variables [Nelson,1989]. The perturbation variables are the 
instantaneous changes from the reference conditions of translational velocities, angular 
velocities, control deflections, and their derivatives. For example, the term 6X/du is the 
stability derivative defining the change in X force with respect to the change in forward 


speed. This derivative can be expressed in terms of a non-dimensional coefficient Cy, as 





follows: 
ox 1 
= ay ——. pal) 
= Somos (2.1) 
where 
dCy 
. ry === (22) 
Cx 6(u/u_) 


oO 


and Q is the dynamic pressure 1/2pV‘. 


Figure 2.8 lists the non-dimensional coefficients used in this model. these coefficients 
are generally broken down into three categories, lateral, longitudinal and control. The 
longitudinal coefficients represent forces effecting the longitudinal axes of the aircraft, 
whule the lateral coefficients represent forces affecting the lateral axes of the aircraft. Non- 
dimensional coefficients, generated in actual aircraft testing, are avaulable for most aircraft. 
By using these coefficients in combination with the dynamics equations, it is possible to 


build a general use flight simulator. 
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Longitudinal Coefficients 


Reference Lift at zero angle of attack 
Reference Drag at zero angle of attack 
Lift curve slope 

Drag curve slope 

Pitch moment 

Pitch moment due to angle of attack 
Lift due to pitch rate 

Pitch moment due to pitch rate 

Lift due to angle of attack rate 

Pitch moment due to angle of attack rate 


Lateral Coefficients 


Cyp Side force due to sideslip 


Cip Dihedral effect 

on Roll damping 

on Roll due to yaw rate 

Cup Weather cocking stability 
Cup Rudder adverse yaw 

Gur Yaw dampening 


Control Coe fficients 


C 
c 


Lbe Lift due to elevator 


DSe Drag due to elevator 
Pitch due to elevator 
Roll due to avleron 
Roll due to rudder 
Yaw due to aileron 


Yaw due to rudder 





Figure 2.8: Aircraft Specification Notation 


1. Aerodynamic Forces and Moments 


Using the non-dimensional coefficients, lift, drag and sideforce are calculated as 


follows: 








2 
cee C ae (Vt+ AVE) pvtvs (2.3) 
L'= pot Cg CqO5he $C, cits + Cy 5,56) SEEAVE e ; | 
2 
(Vt+AVE) pVvts 


ib S (2.5) 





Once lift, drag and sideforce are calculated, these forces are translated into forces 
along the aircraft X, Y, and Z axes as shown in Equations 2.6 through 2.8. The terms Fy , 


F,, and F,7 represent the resultant aerodynamic forces. 


F,7 = —L'cosa — Dsine (2.6) 
F,x = L'sina — Dcoso.—SpsinB (2.7) 
Fay = SpcosP (2.8) 


The aerodynamic moments represent the torque forces about the center of the 


aircraft and are determined in the following equations 


b b pV Sb 
Ly = CupB+ CLpP aye + CL Raw + Cy 5,68 + Cy g,5¢|— as (2.9) 


Gita SV Ee) pVt'Sc (2.10) 
Ma, = Cyto + Catt + CqQsee + Oe ee are + Cose8e oe 


2 
b b pVtSb 2k 
N= /c np + CupP ae +CypR ar Cae Ba + Cy5,51| 5 ( ) 





The forces and moments that result from the above calculations are added to other 


forces and moments at this time (Equations 2.12 - 2.17). 


ee Ax TP Thrust (2.12) 
Fy = Fay (2.13) 
Ee = Fie (2.14) 
eet reas (2.15) 
M = Ma +Mepruset Mayro (2.16) 
Ne neta hace onic (2.17) 
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Engine forces such thrust, torque and gyroscopic effect as well as environmental 
forces such as wind shear can have anywhere from a minor to significant effect on the 
forces and moments along all axes of the aircraft [Roskam, 1979]. However, in order to 
limit the complexity of the model, some simplifications are made. Engine thrust is limited 
to the X-axis only and no calculations are made for torque or gyroscopic effect. This can 


easily be included at a future date if desired and when these values are available. 
2. Calculating Linear and Angular Acceleration 


The total force equations are used to determine the linear acceleration of the 


aircraft [Nelson,1989]}: 


| eee | 

U = VR-WQ- gsin@ + _* (2.18) 
. Fy 2.419 
V = WP-UR + gsinocos® + = (2.39) 
_ Fy 

W = UQ~ VP + gcos#cosO + —— (2.20) 


The total moment equations are used to derive the equations for solving for 


angular acceleration: 


e * 


L = lyyP=L 8-1 Oe RO (2.21) 
M = lyyO+ (xy —I57) PR + 1x7 (P —R*) (2.22) 
N = L 2R-IyoP + (yy —lxx) PQ +1,2OR (2.23) 


However, prior to solving for either P or R, an interim step 1s required: 


La ee Oo eerie (2.24) 

N = N= (i ee PO Ieane (2.25) 
Therefore, 

P = (L"Ipg -N'lyz) 7 (xxdoz —la) (2.26) 


Q = (M-~ (xx ~Izz) PR - 1x7 (P?-R°)) /lyy (2.27) 


Rear (Ngee 171,155 =154) (2.28) 
are the equations for angular acceleration. 
3. Calculating Linear and Angular Velocity 


Linear and angular velocities are determined by numerically integrating the 
accelerations. The Trapezoidal Rule, sometimes referred to as the Modified Euler Method, 


is used. The general method for this integration technique ts as follows: 


x 


P= Pyaat (eae? Ja (2.29) 


a 


- 


where 
P,, = new value of P 
P,,.; = previous value of P 
P,, = current rate-of-change of P (obtained from Euler prediction) 


P,,-1 = previous rate-of-change of P 


dt = integration step size 
4. Updating Position 


Because position updates occur in the world coordinate system, a system must 
exist to convert the aircraft’s linear velocities into changes in world position changes. 
Fortunately, this system was proposed years ago by Leonard Euler [Euler,1758]. By 
knowing the current attitude of the aircraft in world coordinates, the linear body velocities 
are easily converted into world rates by applying the following transformation, which 
effectively rotates the [U V W] vector by the Euler angles[Roskum,1980]. The order of 


transformation is important when utilizing this method and proceeds as follows: 


Us, fi 0 0 lly 


Vol = [0 cosd —sing|| V (2.30) 
Ws 0 sind cosd||W 

Use] {cose 0 sine|| Vo 

W 46 —sin8 0 cos@ Wy 

Uw cosy —siny 0 Voge > 32 
Vw] = {sinw cosy 01} Voe aoe 
Wy Calle 


Integrating the resultant velocity vector, now in world coordinates, by the time 


step of the program, a position update is obtained (Equation 2.33). 


2 x 
eee elena 

_ 33 
Zw} Zw | [Ww 


How Euler angles are obtained is addressed in the next chapter. 


IU. METHODS OF ORIENTATION 


A. INTRODUCTION 


The aerodynamic model generates rotational velocities relative to the fixed aucraft 
body coordinate system. But, just as position updates could only be determined after 
converting linear velocities into world coordinates, orientation updates require conversion 
of angular velocities in a similar manner. Three methods exist for defining the conversion 
of angular velocities to orientations in world coordinates, each having its own particular 


advantages and disadvantages. 


The most popular of these three methods is known as the Euler Method. Using a 
sequence of three angles, the Euler Method provides an intuitive description of aircraft 
attitude in world space [Rolfe,1986]. These angles consist of the familiar azimuth angle y, 
the pitch angle 8, and the bank angle @. The next method, which has become popular in 
recent years, 1s the Quaternion Method. Based on the unit sphere, the Quaternion Method 
provides an elegant method of defining rotations through the use of four parameters. Three 
of the coordinates describe the axis of rotation while the fourth is determined by the angle 
through which the rotation occurs [Shoemake, 1985]. The third method of defining 
orientation is the Direction Cosine Matrix. The direction cosines relate the aircraft body 
axis frame to the world reference frame. Direction cosines, as used in flight simulation, can 
be determined from either Euler angles or quaternions, are utilized for transformations 
between axes, and, alternatively can be directly updated by incremental rotation matrices 


[Paul, 1981]. 


Each method has its own particular advantages and disadvantages and their use 
depends on the application and the implementation. Because NPSNET is a networked 


simulator, the orientation model used must not only render orientations in the world of the 
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aircraft being piloted, but also of other aircraft in the world either flying autonomously or 


piloted remotely across the network. 


B. SURVEY OF METHODS 
1. Euler Method 


The most common method of defining an aircraft’s orientation in world space 1s 
by the Euler attitude angles. Starting with the aircraft’s axis origin aligned with the world’s 
axis origin, the Euler angles specify three successive rotations to bring the world 
coordinates into alignment with the aircraft. The fact that there exist 24 possible ways to 
define rotations. each with potentially different results, means that the order of these 
rotations is important. Euler chose the convention of rotating first about the z axis, then 
about the new x axis and finally about the newest z axis. This convention exists in celestial 
mechanics, applied mechanics, and molecular and solid state physics. The convention used 
in quantum mechanics, nuclear physics, and particle physics, chooses to rotate first about 
the z, then the new y, and finally the new z [Burchfiel,1990]. The convention most often 
used in graphics is standard to aerospace engineers and has been proposed for use by DIS 
[UCF/AIST 1990]. Using the right hand rule, rotations are made, first, counterclockwise 
about the z axis by the angle y, then counterclockwise about the new y axis by angle 9, and 


finally counterclockwise about the newest x axis by angle 6 (Figure 3.1). 
The range of values the attitude angles can take are: 


Ww = £7 OFS 


The aerodynamic model generates velocities in body coordinates. As seen in 
Equations 2.30-2.33, the linear body rates were transformed into world rates by application 


of the Euler angles. What follows is a method for obtaining these angles. 
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rotate in y around Z rotate in 6 around Y 


rotate in d around X 





Figure 3.1: Euler Attitude Angle Rotation 


There is a direct relationship between Euler attitude angles and the angular 
velocity of the aircraft around its body axes [Nelson,I‘ °9]. From this relationship 


(Equations 3.1-3.3), the rates of change of the attitude angle: can be derived 


b = P+ Qsinotand + Rcosotan0 


(3.)) 
6 = Qcos—Rsino (3.2) 
Wy = QsindsecO + RcosdsecO (3.3) 
The inverse of the above equations are 
P = 6—wWsinO (3.4) 
Q = Ocosd + Wsindcosé (55) 
R = —Osind + woos cos (3.6) 


Equations 3.1 through 3.3 are also known as the gimbal equations and are quite 
commonly used in simulation. However, a problem exists when pitch, 8, goes through the 
vertical. That is, where pitch becomes +/-(/2). At that point 6 and w become undefined. 


Implementing a flight dynamics model capable of complete vertical maneuvering. 


necessitates “fixing” the code so a division by zero doesn’t occur. 
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2. Direction Cosines 


In the case of a flight simulation, transforming between body coordinates and 
world coordinates is done quite frequently. A convenient way to represent the 
transformation between two coordinate systems is with the direction cosines. Using matrix 
notation and the direction cosines (a, b, c), the transformation from body to world axes is 


expressed by: 


Aw a, by cil fx 
Yw} = [a by cy] /Y (3.7) 
a a, b, Cc, Z 


XxX, Y and Z represent vectors of any kind, such as force, velocity and acceleration. 
The inverse relatronship, converting world coordinates to body coordinates is the 
transpose: 
a, a, a3) | Xp 


x 
y| = |b, by bY (3.8) 
Z 


Cy Cy €3/| Ly 


In terms of the Euler attitude angles, the direction cosines for the above 


transformations are shown in Figure (3.2) 


a, = cosOcosy b, = singdsinOcos y — cosdsiny C, = cososinOcosy — sindsiny 


a, = cosOsiny b, = sindsinOsiny + cosdcos y Cy = cososin@sin y — singdcos y 


a, = —sin® b, = sindcos® C, = cosdcos9 





Figure 3.2: Direction Cosines in Terms of Euler Angles 
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It should be noted that as there are 12 ways in which Euler angles can be defined, 
there also exist more ways to define the direction cosines. They depend on the coordinate 


system utilized. 


The direction cosines are absolutely necessary for transformations between 
coordinate systems, whether Euler angles or quaternions are used. Direction cosines were 
already used for transforming linear velocities in body coordinates to world coordinates 
(Equations 2.30-2.32). By multiplying those transformation matrices into one matrix, the 
result would be identical to the transformation matrix of Equation 3.7. By using direction 


cosines, the need for determining the intermediate velocities is eliminated. 
3. Quaternion Method 


An alternate method that has gained popularity in the graphics community in the 
mid 1980's is through the use of the unit quaternion. Not a new method, quaternions have 
been around for over acentury. Augmenting the “Four-Parameter Method”, they have been 
useful to aerodynamic engineers for some time and are still the method of choice for 
describing spacecraft orientation [Mitchell, 1965]. Discovered by Sir William Rowan 
Hamilton in 1843 as a result of a search for a successor to complex numbers, quaternions 


provide an efficient means for updating orientations [Shoemake, 1985]. 


a. Defining Quaternions 


There are numerous ways to interpret the quaternion mathematically. The can 
be described as an algebraic quantity, 


w+ixtjy +kz (3.9) 


as a point in three dimensional projective space (w, x, y, Z), aS a linear transformation of 


four space (Matrix), or as a scalar plus 3-vector: 


(w, ¥) ¥ = ik+jy+kz (3.10) 
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The best notation depends on their intended use. The most intuitive approach 
is to view the quaternion as a scalar plus 3-vector (Equation 3.10). However, for algebraic 


manipulation, Equation 3.9, generally becomes more useful. 


A common way of defining quaternion orientation is in combination with 
Euler’s Theorem which states that the orientation of a rigid body can be described as a 
rotation about an axis ¥ by rotation angle ® (Figure 3.3)[Goldstien, 1980]. Constraining the 


vector part to be of unit magnitude, the quaternion becomes: 


©. © 
cos >. Sins (3.14) 


This representation is always of unit magnitude such that: 


wetx ty 427 = | (3.42) 





Figure 3.3: Quaternion Orientation 
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Prior to defining how to rotate a rigid body using the unit quaternion, it is 
necessary to review some of the mathematics associated with the quaternion. For the 
purpose of a flight simulation, an understanding of the multiplication and the quaternion 
derivative is necessary. More comprehensive reviews are obtainable from 


{Shoemake,1985; Golstein, 1980]. 


b. Quaternion Multiplication 


Almost every application involving quaternions makes use of the 
mathematics associated with their multiplication. Sumuar to the algebra associated with 
imaginary numbers, quaternions have three imaginary units, 1, j, and k and are non- 


commutative under multiplication with 


— ie eo = SI (3.13) 
where 
ifok=-ji jk =i = -kj ki = j = —ik (3.14) 


In algebraic notation, the product of quaternion Q multiplied by quaternion Q, 1s 


OOS Wyatt yy + kz) (Ww, Fix, + jy, +kz,) 
=(Ww, —XX,—-Yyy,—2Z,) 
+i(xw, + wx, -—zy,+yz,) (3.15) 
SAMY atx WY) + XZ.) 


+ KC ZW YX) — XY) + wz,) 


and is easily implemented in computer code. Vector notation is more intuitive to understand 


but not as simple to code: 
QQ, = (w.¥) (w,.¥)) = ww, -¥-¥, Ww, tw pV tx, (3.16) 


The result of the above multiplication is a rotation from the orientation represented by Q to 


the new cumulative orientation of Q and Q, in quaternion terms. 
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Multiplication provides a method of orientation extrapolation that can be of 
benefit in anetworked simulation. If Q, represents a finite rotation based on an integral time 
step, and Q represents the cumulative rotation, then a repeated multiplying of Q by Q, will 


result in a smooth rotation across a series of update frames [Burchfiel,1990]. 


c. Quaternions in Terms of Direction Cosines 


One frame of axes (body coordinates) can be brought into coincidence with a 
reference frame by a single rotation D about a fixed axis making angles A, B, and C with a 
second reference frame (world coordinates). The four parameters A, B, C, and D, therefore, 
define the orientation of the aircraft body in world coordinates [Rolfe,1986]. The 
transformation matrix relating the body to world coordinates using these four parameters is 


in Figure 3.4. 


1 | 
DeosAtcosicme =D eet ees Gone 
21 : 2 2 
ae 


1- 2sin7Asin 1 1 l l 
“ ea Bo cosh cin D cos Ey 


mt 1 
2cosAcosBsin =D 2cosBcosCsin? =D 


] 1 
]- Jone Sein’ & B 
2 2 


~ 


l l 
a! I 
ee eer se cos Eine 


l 
2 Gee Weve Cee D 2cosB cosCsin” =D 
2 


| 
1- 2sin*Csin®>D 


l l 
CRU ALD eee +2cosAsin=Dcos =D 
2 2 2 ie 





Figure 3.4: Four Parameter Method 


By substituting the following quatemnion parameters, this method is further simplified: 


] 
do = cos; D 
ell , 
q, = cosAsin =D (3.17) 
| 
q, = cosB sin 5D 
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el 
q, = cosC sin , D 


and the transformation matrix becomes, 


pee 

Xw Go +41 — 92-93 * (4442-993) 2 (4942 + 4,95) | Fy 

Aas Dude? 

Vw] = 2 (4192 +4093) 4079) +492— 43 2 (4293-041) || ¥ (3.18) 
eee? 

: 2 (4443-4092)  2 (4293 +4091) Go - 44 - 93 + 93 


which represent the transform based on the unit quaternion. This result is used for 


converting determining position updates as well as orientation updates. 


d. Derivative of the Unit Quaternion 


To update the quaternion from angular accelerations, the following equations 


are used: 


; J 

do = 75 (q,;P+q,Q+q,R) 

; J 

a 5 (qoP + q,R-q,Q) (3.19) 
; i 

o> = lo eae 

, ] 

ate 5 (qoR + q,Q-q>P) 


Because of the constraint that the quaternion be of unit value and assuming an integration 


step size of less than 1, the above set of equations become: 


, ] 

do = — 5 (4yP +. q,Q + q3R) + Adqg 

, J 

es (qgP + q,R-q,Q) +Aq, (3.20) 
J 

Soames (qyQ + q,P—q,R) + Aq, 


. 1 
G3 = 5 (doR +q,Q-q,P) +Aq, 


where 


X= 1- (qg+q,+93+43) (S220) 
e. Euler Angles from Quaternions 


Many of the auxiliary computations involved with a flight sumulation require 
the use of Euler angles. To obtain these angles from the transformation matrix (Figure 3.2), 


the following method is utilized. 


Because pitch is limited to tn/2, the cos(@) is always positive. As a result, 


obtaining these angles is relatively simple. The elevation angle is derived as follows: 


sin@ = —(a,) (3.21) 


C= asin (a3) 


To obtain the azimuth (y) angle it must be noted that, since the cos(@) is 
always positive, the sign value of a, always reflects the sign value of the sin(y): 
a, = cosOsiny (3:22) 
Therefore: 


a 
“q) * (sign [a2]) Ca 


y = BOOS are 


The roll (9) angle is sunilarly obtained: 


C3 
cos8 





) - (sign [b]}) (3.24) 


6 = acos ( 


C. ADVANTAGES AND DISADVANTAGES 


Quaternions and Euler angles each have their own advantages and disadvantages. 
Euler angles use only three components instead of four to represent orientation. If one were 
to send quatemions over a network in place of Euler angles, as has been proposed for DIS, 


network traffic would increase. However, in cases where angular rates remain constant for 
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long periods of time, by extrapolating orientation updates with quatermions, it would be 


necessary to send an update only when an angle rate change occurs [Burchfiel,1990]. 


Quaternions can be computed directly from the dynamic equations, bypassing the 
computation of transcendental functions necessary in computing Euler angles. If each 
transcendental function costs approximately 20 arithmetic operations then the net cost for 
deriving a rotation update using the Euler Method is 94 operations [Burchfiel, 1990]. This 
can be compared to 42 operations using the quatemion method. In flight simulation, Euler 
angles are necessary for use 1n other simulator functions. This means that, approxunately, 
another 64 operations are necessary, making the Euler Method slightly more 
computationally efficient. However, in many cases it is not necessary to compute the angles 
at each orientation update. Network updates should only occur when a rate change occurs. 
Additionally, radars. gyros and attitude indicators can be updated at slower rates. The 
bottom line is that orientations can be achieved very efficiently utilizing quaternions, 
without calculating Euler angles. When Euler angles are needed for other aircraft functions 
then quaternions become less efficient, depending on how often they need to be computed 


(Figure 3.5). 


OPERATION # EULER OPS # QUATERNION OPS 


CREATING ROTATION 











EULER ANGLE 
CONVERSION 


Figure 3.5: Efficiency Comparison of Euler and Quaternion Methods 
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The most significant advantage of quaternions is that no singularity exists when pitch 


angle (8) passes through +n/2. In the Euler Method, P and R both become undefined due 
to division by zero. Techniques exist, however, for working around this singularity. 
Truncating values as +m/2 is approached will avoid this problem. If the pitch angle is 
truncated at values of 89.99 and 90.01 then a 0.02 degree rotation skip results 
[Goldiez,1991]. Depending on the speed of the program and the rotation rates desired, this 
nay not be noticeable. However, in more fine grained simulations, where a slow vertical 


maneuver is executed, it 1s a factor. 
D. PROPOSED ORIENTATION MODEL 


The use of quaternions 1s desirable for many reasons. Since a dynamics model 1s used 
to drive the piloted aircraft, and the graphics platform utilizes transformation matrices for 
rendering, using the quaternion method offers an efficient method of converting body rates 


to orientations 1n world space. 


Numerous aircraft operate autonomously within the prototype simulation that change 
very little in angular velocity. When this simulator is eventually networked to other 
workstations, quaternions will provide a way to foward interpolate the resulting changes 
in orfentation on the remote computers, thereby eliminating the need for the continued 
transmission of update packets. If updates are eventually needed, the quatemion rates can 


quickly and easily be converted into Euler angles for transmission. 


As the number of aircraft increase in the simulated world, the number of calculations 
necessary for orientation updates begins to multiply. Since only the currently piloted 
aircraft makes use of the Euler angles for additional simulation functions, calculating Euler 
angles is not necessary for all other aircraft in the simulation. Therefore, it becomes more 
efficient to utilize quaternions for defining these rotations, saving approximately 42 


arithmetic operations per update per aircraft. 
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In the piloted aircraft, at the expense of 22 calculations per orientation update, a 
smooth rotation exhibiting no singularity about its axes during vertical maneuvers is 
possible. Future versions of this model will incorporate slow speed, fine grain functionality, 
so smooth changes 1n orientation about all axes is umportant and are best obtained using the 


quatemion. 


pay) 


IV. IMPLEMENTATION ISSUES 


Aside from converting the mathematical model directly to code, other considerations 
must be addressed in the simulation. First, 1t must be modular enough to be easily 
incorporated into the existing NPSNET networked battlefield simulation trainer. NPSNET 
requires that it be operable from a Silicon Graphics workstation using a spaceball 
peripheral for flight control input, and provide update information for transmission to other 
workstations utilizing existing network protocols. Second, it must be able to simulate many 
different types of aircraft and allow interaction with a potentially large number other 
aircraft, either puloted from other workstations or operating autonomously. Third, for future 
use by aeronautical engineers, the system must be designed so that the flight characteristics, 
of whichever aircraft being flown, are based on an aircraft’s dimensional characteristics 


and stability coefficients. 


A. OVERALL SYSTEM LAYOUT 


To satisfy the prototype’s basic requirements, the overall structure of the system 


is designed as shown 1n (Figure 4.1). An aircraft data file makes it sumple to create new 


Aerodynamic 
Model 
(Piloted Acft) 


Euler Angle 
Calculation 


Flight 
Data 
Structure 


Onentation 
Model 


Aucraft Data 
File 


Non-Dynamic 
Model Display 
(Autonomous Pipeline 
Aircraft) 





Figure 4.1: Prototype Flight Simulator Basic Structure 


aircraft, position these aircraft, and designate their handling characteristics. It is also be 
used to initialize the flight parameters of the autonomous aircraft. The flight data structure 


contains information describing the current state of the aircraft and its design 


30 


specifications. The aerodynamic model is used for updating the piloted aircraft’s body 
rates. The non-dynamic model also outputs a set of velocities, but unlike the dynamic 
model, it determines these velocities in accordance with a predetermined script designated 
by the implemented. The orientation model converts body rates to positron and orientation 
in world space. Euler angles are then determined for the puloted aircraft and are needed to 
update cockpit displays. Autonomous aircraft do not require the calculation of Euler angles 


and bypass this function. 
B. THE AIRCRAFT DATA INPUT FILE 


Data records within the input file are divided into two categories, flight records 
and aircraft specification records. The flight record contains information describing the 
position of an aircraft and its initial flight parameters (Figure 1.2). The specification record 
contains the dimensional characteristics and stability coeff. ients describing a particular 
aircraft. By manipulating the data in the specification record, one can change an aircraft's 
basic design as desired. To link an aircraft to a particular set of specification data, all that 
is necessary 1s to match the “type aircraft” information within the two records. This system 


allows one specification record to be used for an unlimited number of flight records. 
C. THE FLIGHT DATA STRUCTURE 


The flight data structure provides a global source of information on the state of each 
aircraft in the simulation. The information contained here is necessary for operation of the 
aerodynamic and orientation models, and for updating cockpit displays (Figure 4.3). Note 
that the fourth item in this structure is a pointer to another data structure containing aircraft 
specification data. Not shown, this structure contains information identical to that found in 
the specification record of the data input file (Figure 4.2). Maintaining flight information 
in two separate files saves some storage space by allowing more than one aircraft use a 


single set of specification data. 


Flight Record 


| 

I 
200.0 
-950.0 
300.0 
0.0 
0.0 
0.0 
a0) 

| 


Specification Record 


4 

I 

2 
260.0 
10.8 
546.0 
8090.0 
25900.0 
29200.0 
1300.0 
0000.0 
8000.0 
0.03 0.3 


/id number/ 

/type aircraft/ 

/airspeed/ 

/posx/ 

/alutude/ 

/posz/ 

/heading/ 

Anitial angle of bank/ 

Ainitial gforce/ 

/status:0-piloted,1-levelturn,2-g controlled turn/ 


/type aircraft-- A4/ 

/jet or prop jet:1 prop:0/ 
/o/ 

[S/ 

/c/ 

/m/ 

/\x/ 

[ly/ 

/\z/ 

[xz/ 

/max thrust/ 

/mil thrust or horsepower/ 
/CDo /CDa 





0.28 3.4500 0.72 0.36 /CLo /CLa /CLq /CLda /CLile 

0.0 -3.6 -0.38 -1.1 -0.5 /CMo /CMg /CMa /CMda /CMde 

-0.98 0.17 /CYb /CYdr 

-0.12 -0.26 0.14 0.08 -0.105 = /CLb /CLp /CLr /CLda /CLdr 

0.25 0.022 -0.35 0.06 0.032  /CNb/CNp /CNr /CNda /CNdr 

0.2618 -0.5236 0.5236 /deflection limits of rud, ail, elevator (radians)/ 


Figure 4.2: Data Input Records 
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int ID: --number assignment 

int type: --aircraft type 

int Status; --piloted:0 or autonomous level turn: 1 climbing tum:2 
typeptr Tptr; --pointer to aircraft specification data structure 
float Forces|[ 3}: --forces in X,Y,Z dir 

float Torques{ 3]; --torques around X,Y,Z axis 

float linear_vel[3] : --velocity in X,Y,Z direction 

float angular_vel{3]; --angular velocity around X,Y,Z axes 

float linear_accel[3]: --linear accelerations 

float angular_accel{3]: --angular accelerations 

float sideslip: --Sideslip or beta angle 

float ang_atk: --angle of attack 

float d_ang_atk; --angle of attack rate 

float Lift: --total lift 

float drag: --total drag 

double Q{4}: --quatemion 

Matrix at --direction cosine matrix 

float euler_angles[3]: --euler angles angles in radians - yaw, pitch.roll 
float pos{3]: --world position in X, Y, Z 

float ref[3]: --look direction 

float vel_world{3]; --velocities in world position - X, Y,Z 

float gfor: --amount of g force 

float rpm: --engine mm 

float elev: --flight control positions 

float eltrm; --elevator trim 

float --rudder position 

float aul: --aileron position 

float --throttle position 

int --flap position 

int --landing gear position 





Figure 4.3: Aircraft Flight Data Structure 


D. FLIGHT CONTROL INPUT 
1. Throttle 


The throttle un an aircraft is the pilot's primary means to control] the engine. As 
such, a mapping of throttle position to engine rpm must be devised that incorporates delays 
associated with engine spool-up characteristics. In Jarge, high speed simulators, rpm and 


engine data is retrieved from engine-specific lookup tables. Because tables such as these 


are engine specific, a sunpler, generic method was devised. By keeping track of the tast 
throttle position, the following formula is utilized: 

Arpm = (T,a—- They) Pat (4.1) 
where 

T = throttle position 

dt = delta time 

¢ = engine spool-up delay factor 

Since applications using this simulator include both jet and propeller aircraft, it 

was decided that allowances should be made for the differing characteristics of the two 
types of engines. Jet engines are generally rated in terms of thrust (lbs), while propellers are 
rated in terms of horsepower (ft-Ibs/s) [Anderson, 1989]. Since the aerodynamic model 
uses thrust in terms of Ibs, it can use the data for jets durectly. However, propeller driven 
aircraft require the following conversion: 


550nHP p 


oe? (4.2) 


where 
n = propeller efficiency (usually around .8) 


HP = engine rated horsepower 


Pe. density altitude ratio where p_ 1s the density at sea level 
p 9° 


0 


2. Ajileron, Elevator and Rudder Control 


A nonmal aircraft stick exhibits two degrees of freedom, left-right for aileron 
control and back-forward for elevator control. Therefore, control inputs from the spaceball 
were limited to these directions. The maximum deflection of the control stick is information 
entered via the specification records. It is a simple procedure to read deflection data from 
the spaceball and linearly map it to a control deflection somewhere between +/- max 


obtainable deflection. 
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Rudder input is handled in a simular manner to aileron and elevator control. 
Maximum rudder deflection 1s entered via the specification records so that its position can 
be linearly determined from the amount of rudder control input. Rudder input 1s via the 
keyboard using the left and right arrow keys. Although not the same as using rudder pedals, 
separating the rudder input from the contro] stick inputs forces the operator to be concerned 


with the issue of rudder and aileron coordination. 
3. CONTROL OF NON-PILOTED AIRCRAFT 


In the current system, an autonomous aircraft can exhibit one of two behaviors, 
depending on the value of the “status” data in the flight data record. The first behavior is 
that of a level turn. The initial roll angle, also a part of the flight data record (degrees) 
provides the seed for calculating the desired body rates. From the bank angle, the amount 


of g force necessary to maintain level flight at this angle of bank 1s [Anderson, 1989]: 


| 
G-force = cosé (4.3) 


from the g-force, the rate of turn (ROT) in degrees/sec can be determined: 


g/G-force ang 


Se (4.4) 


ic 


where g is gravity (32.2 ft/sec’). A level tum requires a balance between the angular 
velocities Q, pitch rate, and R, yaw rate. Therefore, from Equations 3.4-3.6: 


Q = ROTsino 


(4.5) 
R = ROT cos 


The second behavior is that of a non-level turn. This turn is based on the initial 
angle of bank and g force entered in the data input file. The yaw rate, R, is calculated the 
same as for a level turn. The pitch rate, Q, is also computed using these formulas but, 
instead of using the computed g force for a level tum, it is calculated using the amount of 


desired g force entered via the data file. The calculations necessary to compute the angular 


aS 


velocities are completed during program initialization and do not impact on the run tume of 


the simulation. 
E. SPEED OF AERODYNAMIC MODEL 


The aerodynamic model exhibits a tendency to “blow up” if the time step between 
calculations becomes too great. This becomes evident when the aircraft speed and rotation 
rates increase. The solution to this problem is to run the aerodynamic model at a faster rate 
than the rest of the system. The trick is to determine how fast. “Smart” integration schemes 
can be found in the literature that increase or decrease the number of time steps according 


to the amount of numerical divergence present in the integrated values [Press.et al, 1990]. 


A simple method 1s used in this program to solve for this problem. Because there exists 
a direct relationship between an aircraft’s speed and its tendency to produce a numerical 
instability, the time step was adjusted in direct relationship with the aircraft’s airspeed. 
Prior to updating body rates in the aerodynamic model, aircraft airspeed is measured and a 
the number of time steps is calculated (Figure 4.4). The timestep factor used in the current 


program 1s 0.03 seconds. 


aerodynamic model 
read aircraft velocity 
time step for aero model = timestep_factor * velocity 
for computed_time_step loop 


do aero calculations 
update aircraft state variables 
end loop 
exit aerodynamic model 





Figure 4.4: Algorithm for Computing Time Step 


Increasing the time step does not decrease performance of the overall system. The 
speed of the aerodynamic model in the current unplementation ranges from 17.6 ms to 18.8 
ms for time steps ranging from two to 160. Running at approximately 120 ms, the graphics 


pipeline definitely remains the limiting factor in this system. 
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V. CONCLUSIONS AND FUTURE WORK 


The techniques presented in this paper have proven to be an effective method of 
implementing a graphical dynamic flight simulation on a matrix-based graphics computer 
in real-time. Like most research and academic projects, this aircraft sumulator is structured 
to allow for the addition of more detailed functionality. Current work includes the 


integration of a weapons delivery system and avionics suite. 


The orrentation model functions developed during the course of this paper will] become 
part of the C program library at the school, thereby providing an alternate and more flexible 
tool for manipulating solid objects tn a graphical environment. Integration of this 


orientation model into other dynamic simulation systems is « .o under investigation. 


APPENDIX A 


[RRKKS KHKKK KKK KK KKK KEK KKK KEK KKK KKK KKK KE KKK KE KKK KKK KKK KEKK KK KK KK KKK KEKKKKKK KKK KKK / 


/* FILE: aero.c * / 
/* AUTHOR: Joe Cooke oe! 
/* DATE: 1/9/92 x / 


/* DESCRIPTION: This file contains the mathematical model for all */ 


/* aerodynamic calculations. * / 
/* Input: aircraft state structure (P) * / 
/* aircraft specification stucture (T) * / 
/* Output:updated state structure (P) - 


[ORO III III III IOI III III IGIGI III III III IOI IOI III IOI IOI IO III IIIT / 
aero model (acftptr P,typeptr T,float deltatime) 


{ 

float 
float 
float 
float 
float 


0, OS, 0se, Geb, e200, 220, 4,4 mn 
totave l: 
densaltratio,densalt; 

aoa, Cos aoa, sin aoa; 
cos_sideslip, Sin_sideslip; 


/*misec. variables*/ 
/*total velocity*/ 
/*atmosperic data*/ 
/xangle of attack data*/ 
/*sideslip data*/ 


float For thrust x: /*thrust force*/ 

float ailpos,rudpos, elevpos; /*control positions*/ 
float Itemp, Ltemp, Ntemp; /*temporary terms*/ 

float du, dv, “dw, dp; aq, ar. /*newly computed accels*/ 


float 
float 
float 
float 
Float 
Eikoat 
float 


diff, 

loopdelta; 
vworld[3], vwer [sti 
Valrerare| 2), 
hpower; 

max thrust; 
loop_iterations; 


Tne 2st oe ee 
Matrix A MATRIX? 


[RKKKKKEKKKKKKKKKK EASE 
it (e 
PrLUAtt (“****error: 


== NULL || == NULL) { 


return; 


} 


/*temp var for rpm determination*/ 
/*deltatime / loop iterations*/ 
/*vel in world inertial frame*/ 
/*vel in aicraft frame*/ 
/*horsepower*/ 

/*thrust*/ 

/*number of aero calc iterations*/ 
/*index and loop vars*/ 
/*temporary matrix* 


for valid SELCUCCULESKKAKKKKKKKKKKKKKK KK KKK / 


null pointer Sent to Aero.c****\nvn 


/****some initial assignments****/ 


no =| Sm. 
g = GRAVITY; 
w=m * GRAVITY; 


P->gfor = P=>lift/w; /*calculate current G’s*/ 
if (P-—>u < 10.0) P=>u = 10.0; 
tot vel =» feqrt (P—>v*P=>y et Po>w bow ete Pe), 


[xkkkkakkKKKKKKatmosheric density CalculationsKr KKK KKK KK KKK RRR K / 
if (P=>posy >= 0.0 46. P—=posy < MAX TAIT) 

alt = (int) (P->posy/ (ALT_INCREMENT*FT TO METERS) ); 

densalt = densaltarray([alt]; 

densaltratio = densalt/DENSITYSL; 

} 


else{ 
densalt = DENSITYSL; 
densaltratio = densalt/DENSITYSL; 
} 


[KKKKKKKKKKKKKKKCAaAlCUlate new CPM BRK KKKKKKKEKKKKKKKKKKKKEKEKKKEK K / 
diff = (P->thro -— P=->rpm) /(2.0/deltatime) ; 
P=<>rpm = P—>rpm + diff; 


[OOO IK Calculate new thrust * RII IRI IO IIR I IR / 
if (T->class == 0){ /*if a prop plane*/ 

hpower = PROP_EFFIC * T->horsepwr; 

max thrust = (hpower * HPR2THRUST)/ tot vel; 
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Pome noiacek = mas eenruss ~ P—>rpm * densaltratio; 


else{ /*if a jet plane*/ 
ee Ae Conor eh hrust x o=—(T—->max*(P—>rpm)) * densaltratio; 
Sew noretinniget xe (lm * (P—>rpm) jes densaltratio; 
) 


[xrxxraxkekdetermine control differential gt ** eK KKKKKKKKHKKE RK / 
meepos = {(P—>ail * T->def #il); 

mpacapos- = (P—>rud * T->def rud) ; 

elevpos = (P->elev * T->def_elev) + P->eltrim; 


fRRKKKKEXKKKKadjust time step for aero calculationst****e*rKkx / 
Pompei cerations = (float) ((int) (tot_vel* TENTH OF AIRSPEED FACTOR) ); 
loopdelta = deltatime * 1.0/loop iterations; 


[RKKKKKKKKKK]LOOp to reduce time step for dynamic model*****/ 
Poctest = Ofixt<((int)loop iterations) ;ixt ++) { 
if (P=>u < 10.0) P=>u = 10.0; 
cae Vet agit (P—>v*P=>yv + P=>whP—>w + P—>u*P->u) ; 
P-—>sideslip = fatan((P->v/P->u) ); 
sin sideslip= (sin(P->sideslip) ); 
cos sideslip= (cos (P-—>sideslip) ) 


e 
, 


aoa = fatan( (P—>w/P=—>u));/* + (5.0/RAD TO DEG) */ 
cos aoa = (fcos(aoa)); 

Sin aoa = (fsin(aoa)); 

P->d_ ang atk = aoa -— P->ang atk; 

P->ang_atk = aca; 


OF. densalt*tot vel*tot vel; 
GSa= 30 *T—>S? 
QSc = QS*T=->c; 


OSb = QOS*T—>b; : 
ez tec, (2.0*tot vel); 
b2up =) T—>b/ (2.0*tot vel); 


fur ee calculating aerodynamic LOrceskkRKKKKKKKKKKKKKKKKKR / 
Bite) tt —>CLo + 

tea clkag ro mang atk + 

Resenae bog) * C2 + 

Tarobdalphas* P->d ang atk *'c2u + 

i=-Cude *seleyvpos ) * QOS; 


Fo->drag = (1f->CDo + 
Lace wee eng atk ) * OS; 


P—>sideforce = (T=->CY¥b * P=->sideslip + 
i=rordr * rudpos ) * QS; 


P->Fxé = (P->lift * sin aoa - 

P->drag * cos aoa - = 

P->sideforce * ae cle eta 

Eo -hy.= (P=>sideforce * cos sideslip); 

aoe eee litt * Cos. aoa — P—->drag * sin,saoa); 


[ktKKKKKKKKCalCUlating aerodynamic MOMENtS*KAKKKKKKKKKKKRKE RRR / 
P—>L = (T=->CLb * P->sideslip + 

hoa Clip ~ P—>p = b2u + 

taecbr * Pe>r * b2a + 

T->CLda * ailpos + 

foe chars* rudpos. ) * OSb; 


P->M = (T->CMo + 

T—>CMa * P->ang atk + 

T->CMq * P->q * c2u + 

=7Cida * P=-deang atk * c2u + 


T->CMde * elevpos ) * QSc; 


P=>N | (T=>CNb * P—>sideslip + 
T->CNp * P->p * b2u + 


oe 


T=>CNr * P=>r * bZu 4 
T—>CNda * ailpoes + 
T->CNdr * wudpes )) = Osnp; 


/******determine total moments and forces*®*® KKK KKKKKAKKKKEKKKE KS / 
P->Fx = P-=->Fx + For thrust. 


/**k*k*temporary calculations ® * 677 IORI III III II IOI III 
Ltemp = P->L + T->Ixz * P=>p * P->q = (TSize tee eee 

Ntemp P=>N — (TH=>Iy =— T=>1x) * Poop * Pasa ae ee ee eee 

Itemp (T=>Ix.* T=siz i — f-Fisz “sie is 


/xxxkkkdetermine Linear and Angular Acceleratlions* *KKEKKKKKKEKKRKEKKKKAKKERE / 
du = (P=->v * Pe>r — P=>w * P->q + P->Fx/m - g * P=>sptch ); 
dv = (P=>p * P=>w - P=>r * P=>u + P=>Fy/m + g * P=>sroll * P=>cptch); 
dw = (P=>u * P=>q — P=>v * P—->p + P=>Fz/m fq 8 ClO) le See coe) 
dp = (Ltemp * T—>Iz + Neemp * 1-51.2) 7 tomo, 
dq = (P=>M = 9 ((T=>ix — f-7iZie 2. Pe ee 
- ((P->p * P->p - P=->r * P=>r) * T=>Ixz)) / T->Iy; 
dr = (Ntemp * T->Ix + Ltemp * T->Ixz) / Itemp; 


/***determine velocities with trapizoidal integration****/ 


P—->u = P->u + (( P->udot + du )/2.0) *loopdelta; 
P->v = P=>v + (( P->vdot + dv )/2.0) *loopdelta; 
P->w = P->w + (( P->wdot + dw )/2.0) *loopdelta; 
P—->p = P->p + (( P=>pdot + dp )/2.0) *loopdelta; 
P—>q = P=—>q + (( P=>qdot + dg )/2.0) *loopdelta; 
P—->r = P=>r + (( P->rdot + dr )/2.0) *loopdelta; 


[*k*KK*KKKKreCOrd Last acceleratLom* KH KIKI RIKI II Ht KH / 
P-—>udot du; P=->vdot = dv; P->wdot = dw; 

P=->pdot = dp; P=>qdot = dq;P-—>rdot = dr; 

} /***end loop***/ 


{/**k*x**kincrementally Totate GTUAC|SrnLONnk KKKKAKKKAKKKEKKKKKKKEE / 
increment _quat (P->p, P—>q, P—>r, P->Q, deltatime) ; 


/***k*ypdate rotation matrix from quaterniom: ~="**" <== === 
rotation matrix from quat (P=>90, P=>H); 


[**kk*eextract euler angles From Mat Link KX KKKKKKKKAKAKAKERER / 

euler _ angles from_matrix(&P-—>sroll, &P—>croll, &P-—>sptch, &P-—>cptch, 
&P=>shdg, &P=>chdg,P-—-h) matrix, 

euler _ang frm rot matrix(&P=>r0ll, &6P—>pteh, «P—--ndq, P—>r, 

P=>ptch = P=>ptch * RAD TO DEG; 

P->roll = P=>roll * RAD TO DEG; 

P->hdg = P->hdg * RAD TO DEG; 


/***kcalculate new world velocity ****** x4 ens ee RRR R RRA ee 
vaircraft(0] = P->u; 

vaircraft(1l] = P->v; 

vaircraft (2] = P->w; 

transpose matrix(P—->H,A_ MATRIX) ; 

post_mult_matrix_by vector(A_MATRIX, vaircraft,vworld) ; 


vwor[{0) = vworld[0) * FT TO_METERS; 
vwor(l] = -vworld(2]))* biome ibes: 
vwor[ 2] = vworld(1} * FT TO METERS; 


/*xxkcalculate new world POSACION*** **xKAKEKAARAKARARE AAA 
P->H(3](0] = P->posx += vwor[(0] * deltatime; 

P->H(3](1] = P->posy += vwor(1] * deltatime; 

P=->H(3]J (2) P->posz += vwor[{2) * deltatime; 

P->H(3)] (3] 0 


[*xekkkkkcalculate world reference Point **XKKKHKKHKKK KEKE / 
P=>refx = P—>post + (2->eptch) = Y= end). 

P->refz = P->posz + (P=>cptch * P->shdg); 

P->refy = P->posy + P=->sptch; 
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APPENDIX B 


[RR KKK KKK KKK KKK RK KR KK KK KR KKK KK RK IK KKK KKK KKK IKK KK KK KKK KKK KKK KR KK IK RK / 


/* FILE: quaternions.c * / 
/* RUTHOR: Joe Cooke x / 
/* DATE: 1/9/92 * / 
/* DESCRIPTION: This file contains the functions for the quaternion * / 
* 


/* rotation and matrix operations 
[RIO OO IO OOO IOI III IOI IO IO IO IO IOI IO IO III II I I aK / 


{/***multiply quaternion Qr by Ql and return the product in Qr***/ 
mult quaternion (Q1,Qr) 
double Q1[4],Qr[4]; 


{ 
double n_factor,q; 


fmeeaeeos— 1-07 ((Or[0) *Or[0)])+(Or[1] *Or[(1]) + (Or[(2] *Or[2])+(Or[3] *Or[3])); 


Q1(2]*Or[2] - Q1[3]*Qr[3]) 
ote | onl> le Clls]*Or[2))) 
Cis eatin Cl [1)*Or[3))) 
Oe Orie Cbl2)*Orl1)) 


Peei—1 factor* (O1[0)]*Or[(O] — Q1[(1] *Or[1] 
Seaiyean Lactor* (Ol[i]*Or[0)] + Q1[0]*Or[1) 
eee | factor~ (01 [2] *Or[0] + Q1[0] *Or[2] 
feo p=ne factor (O1[3]*Or[O0] + Q1(0) *Or[3] 
} 


f 
¢ 
? 


+++ 


* 
? 


/x*xconvert roll, pitch and yaw angles to a single quaternion*****/ 
ae wsetiul for intitalizing a quaternion with initial orientation***/ 
euler to _quaternion (roll, pitch, yaw, Q) 

Eroee roll, pitch, yaw; 

double Q[4]; 


Sonblelt,p, V,cosr,cosp, cosy, sinr,sinp, siny; 


r = (( (double) roll/2.0)/RAD TO DEG); 

p = (double) ((pitch/2.0)/RAD TO DEG); 

y = (double) ((yaw/2.0)/RAD TO DEG); 

eoer = Ccos(r); cosp = cos(p); cosy = cos(y); 
Seat = Sin(x)? Sinp = sin(p); siny = sin(y); 
Q(9)] = (cosr*cosp*cosy)+(sinr*sinp*siny) ; 
Q[1] = (sinr*cosp*cosy)—(cosr*sinp*siny) ; 
Q[2]) = (cosr*sinp*cosy)+(sinr*cosp*siny); 
Q[3] =(-sinr*sinp*cosy)+(cosr*cosp*siny) ; 


/xx*increment a quaternion from incremental body angular rates******/ 
/***p=roll q=pitch and r=yaw Celt a grrr eenk nee ke ekkKKKK KKK KK RK KK KR / 
increment quat (p,q,r,Q,deltatime) 

Post p,q, x; 

double Q[4]; 

float deltatime; 


Couple factor, qql,qq2,qq3,qq0,dq0,dql,dq2,dq3,p2,q2, x2; 


ago = ©0[0)]*Q(0]; 
aaee= O[1)*O[1); 

=—O02)*O[2)% 
aaee= ©(3)*O(3); 


pzZ = (double) (p- 0. 5*deltatime) ; 
q2 = (double) (q*0.5*deltatime) ; 
r2 = (double) (r*0.5*deltatime) ; 


Eactor = 1 .0-(qq0tqqlit+qq2+qq3) > 


dqO = -(Q[1)*p2 + Q[2]*q2 + O[3]*r2) + Q[0]*factor; 


4] 


dql = (Q[0]*p2 + O[2]*r2 - Q[3]*q2) + Q[1]*factor; 
(Q[0])*q2 + O[3)]*p2 = Of) *r2Z) + O[2)]“*tacror 


dq3 (OG) *nZ +O 1)" 4q2-— O21 *p2) + Oe tae 


Q. 
Q 
N) 
| 


Q[0)]. = "O]01) + aq0- 
O[1) “= Obi + dal; 
O[2] = O12)) 4 daz; 
O(3) = Oe) 2 4-aq3- 


/***matrix for conversion from worldrates to body rates***/ 
rotation matrix from_quat (Q,M) 

double Q[4]; 

Matrix M:; 

double £,wf,xf,y£, 2f,*xE,yVt,oct, xvi, Xztp vets wy twat, 


f = 2.0/(Q[0] *O(0) +Q[1) *QO[1] +Q[2] *O[2] +Q[3]) *Q[3])): 


xf = £ * Of[1); 

co £ Xo: 

Zf = £ * Of 3)- 

Sf = “PFO 

VVE = yEror 7s 

zeL = zE*OlSl- 

MVE = 2270 (24). 

27 = or Oo 

yzi = VyE*Ols)> 

wxf = ximO(0]: 

wyf = yf*Q[(O); 

wzf = z£*Q(0); 

M(O) [0] = (float r= e4222)); 
M(O)[1] = (float) (wzf+xyf) ; 

M[O} (2) = FGl lca) (zt en ee, : 
M(1] (0) (= tileae) (xy i—-w2i) 

M(1) (1] c= (£loeat) (1.0=(2xF+4+z2z£)): 
M[1) [2] 9m “ileat)y (vzitwxt) 

M[2] (0) = (float) (xzf+wyf); 

M(2) (1) = (float) (yzf-wxf) ; 
MiZ)[2} = (float eo et tyvew 
M[3) [0] = 0O.; M[3] [1] = 0.7 M([S)(2) = 0.) UC eS: 


) 


/***this function “extracts” sine and cosine data from a cosine matrix 
represented in conventional Lorm, KR KKK KKK KKK KKK KK KK KK KK KK KK KK / 
euler angles from _matrix(sroll,croll, sptch, cptch, shdg, chdg,M) 
float *sroll, *croll, *sptch,*cptch,“*shdg,enag, 
Matrix M; 
{ 


Pleaewsp,cp,s0,cn, SV, cy: 


sp = -M[0] [2]; /*sinpitch*/ 
cp = faqrt (1.0-(sp*sp))- j/*ecspitch* 7 
if (cp == 0.0) ( 7*if pitéeh at 47— 90deq*/ 
sy = 0.0; /*sinyaw */ 
cy. = 1.0) /*cosyaw */ 
sr = M[1] [0]: J ®sanroll*/ 
cee M[Muie /*cosroll*/ 
} 
else ( 
Sy eo ny co) /*sinyaw */ 
ey 2e MiO) (O07 (ep) > /*cosyaw */ 
sx = M[(1])[(2]/ (cp): /xsinroll */ 
era MZ) (ep: /xcosroll */ 


} 
ksptch = sp;*cptch = ¢p; *sroll "Segsr?*croll = cr; *shdceye +);  chagu-ac,., 
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/***The following function will return the roll, pitch and yaw from a standard 
transformation matrix.This routine should work for a transformation in a coordinate 
system from body coords where the x y and z axes move in pitch roll and yaw 


respectively, to world coordinates where X Y and Z point North, East and Down 
steric eat me eee a er AAR AR AR RARRAE RRA RERKE RE / 


euler ang _frm_rot _matrix(roll, pitch, yaw, matrix} 
float *roll,*pitch, tyaw; 
Matrix matrix; 


Bloat cOSspitch, Sinpitch, roll fact, yaw fact; 


/**determine pitch angle in radians and cospitch**/ 
Sempi1ctch = —matrix([0) (2); 

Soepetch = fsqrt (1.0 - (sinpitch * sinpitch)}; 
pouten = —(fasin(matrix([0] [2))); 


/xkdetermine roll angle in radians**/ 
Bolle fact = matrix(2) (2) /cospitch; 


if (cospitch != 0.0} 
{ 
Mime mats x1) [2] < 0.0) *roll = —-(facos (roll fact)); 
else 
peel = facosi(roll fact); 
} 
else 


cro! = 0.0 


/**xdetermine yaw angle in radians**/ 
yaw fact = matrix[0]) [0]/cospitch; 


et (cospitch i= 0.0) 
1 
if (matrix[1)[0] < 0.0} *yaw = -(facos (yaw fact) ); 
else 
yewe=sceacos (yaw fact) ; ; 
} 
else 


*vaw = 0.0; 


/***this routine assumes matrix is in graphics convention form. The robotics convention 
requires an additional matrix transposition prior to utilizing this functionk *** kek x / 


euler to acne axis convent (M matrix,W _Mmatrix) 
Matrix M matrix, W matrix; 


{ 


Wematrax([O0) [0] =M matrix[0] [0]; 

W_matrix[1])[0] = M matrix[1] [0]; 

Nematrix(2) [0] =M [| matrix[2] [0]; 

Womatrix[0O][1] = -M matrix[0] [2]; 
Womatrix(1)(1) = -M matrix[1]) [2]; 
Womatrix[(2J][1] = -M_matrix[2] [2]; 
Womatrix(0][(2] = M matrix[0][1]; 

W_matrix(1])[(2] = M matrix[1] [1]; 

femeacrss (2) (2) = M matrix(2) [1]; 

W_matrix[(3]) [0] = M matrix[3] [0]; 

Sameera (S) fl} = M matrix[3] (1)? 

Bemetrix~(3)[(2) = M matrix[3] [2]; 

Perecri (3) (3) = 1.0; 


} 


/**kReturns a vecter result ®** RRR RRO OOOO III IOI TOOT IO IOI / 
void post_mult_matrix_by vector(Matrix M matrix,float * V_array, float * Result) 
{ 

Result[0] = V ee oy, CO) hs matrix[0])[0]) + 

Vv “array [1] *M_ matrix[O][1] + V Arey | 2) *M matrixo 2]; 

Result [1] = V peroay (0) hl matrix Peo 

Pmenray (1 |)*M - matrix[1][1] + V waczay (2) *M mater) (2); 

Result[2] = V ar ray (0) >M matrix[2])[0]) + 

- array (1) *Me matrix[2][1] + V warrvey | 2) Mematrix (2) (2); 
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APPENDIX C 


A. INTEGRATION OF AERODYNAMICS INTO NPSNET 


To complete the integration of the aerodynamic model into NPSNET the following 
steps must be accomplished: 
«Modify the vehicle model structure to fully describe a dynamic aircraft. 
¢Create a method to assign aircraft characteristics on program initialization. 
¢Add cockpit display as the appropriate user interface. 
¢Modify input devices for aircraft control (simulate aircraft control stick). 
¢Add aerodynamic model for control of the puloted aircraft. 


¢Add orientation model based on quaternions for aircraft rotation. 


B. NEW VEHICLE MODEL STRUCTURE 


The current NPSNET has one model for all types of vehicles. Tanks, helicopters, jets 
and jeeps all move under the same sets of rules, and only one basic data structure for all 
vehicles was previously created. Because air vehicles operate in a different dynamic 
environment than ground vehicles, it is necessary to create an additional data structure for 
their description. A logical way to add additional data for an aircraft vehicle is through the 
use of a union data structure. Imbedding a union data structure containing the new 
aerodynamic model within the original data structure allows the insertion of additional data 
needed for a more complex model, without having to re-engineer the entire program with 


a completely new data structure (Figure A3.1). 


... (attached to Old Data Structure For Vehicle Model) 
union { 

DEFVEH defaultveh: 

AIRVEH airveh: 


GNDVEH gndveh: 
SHIPVEH shipveh: 
SIMNETVEH simtype: 
} veh; 





Figure C.1: Union to Be Inserted into Current Vehicle Data Structure 


44 


Two data structures are added to handle the aerodynamic model. Figure A3.2 shows 


the structure containing the vehicle state information and Figure A3.3 shows the structure 


containing the vehicle specifications. To connect the two structures, the pointer, 


TYPEPTR, was created. 


typedef FLIGHTSPECS * TYPEPTR: 


typedef struct 
( 
int totalacft: 
ebeP VR Tptr: 
float torc[ 3]: 
float torg{3]: 
float lin_vel{3}: 
float ang_vel[3]: 
float lin_accel[{3}: 
float ang_accel[3]; 
float sideslip: 
float ang_ atk: 
float d_ang_atk; 
float lift; 
float drag: 
double quat|[4]: 
Matnx h_matmrix: 
float sroll: 
float croll; 
float — sptch; 
float = cptch; 
float shdg: 
float chdg; 
float euler_angles[3}; 
float vel_world[3}; 
float hdg: 
float gfor; 
float rpm; 
float elev; 
float eltrim; 
float rud; 
float ail: 
float thro; 
int flaps; 
int gear; 
} AIRVEH; 





/*track total number of aircraft in simul*/ 
/*pointer to aircraft type information*/ 

fF ibs) Forces in X, YZ direction*/ 

/* ft-lbs: Torques around X,Y,Z axis*/ 
/*ft/sec: linear vel(u,v,w) in X, Y,Z dir*/ 
/*rad/sec: ang vel(p,q.r) around X, Y.Z axes*/ 
/*ftps2: linear accels in X, Y.Z direction*/ 
/*deg/sec2:angaccel(p,q.r) around X,Y,Z axes*/ 
/*sideslip or beta*/ 

/*angle of attack*/ 

/*angle of attack rate*/ 

/*total lift*/ 

/*total drag */ 

/*quaternion*/ 

/*direction cosine matrix of world position */ 
/*sine of roll*/ 

/*cosine of roll*/ 

/*sine of pitch*/ 

/*cosine of pitch*/ 

/*sine of heading*/ 

/*cosine of heading*/ 

/*euler angles in radians - roll,pitch,yaw*/ 
/*meters/sec: velocity in SGI coords X,Y,Z*/ 
/*world heading*/ 

/*g force exerted from back stick*/ 
/*fraction: 0.0 1.0 engine rpm*/ 

/*position: -1.0 to 1.0*/ 

/*elevator trim*/ 

/*position: -1.0 to 1.0*/ 

/*position: -1.0 to 1.0*/ 

/*position: 0.0 to 1.0*/ 

/*boolean: 0-flapsup 1-flapsdown*/ 
/*boolean: 0-gearup 1-geardown*/ 


Figure C.2: Data Stucture Containing Aircraft State Information 
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typedef struct 
{ 
int 
int 
float 
float 
float 
float 
float 
float 
float 
float 
float 
float 
float 
float 
float 
float 
float 
float 
float 
float 
float 
float 
float 
float 
float 
float 
float 
float 
float 
float 
float 
float 
float 
float 
float 
float 
float 
float 
flout 
float 


type: 
class: 
span; 
sarea; 
cord: 
mass; 
Ix: 

ly: 

jeep 

[xz: 
horsepwr; 
niax: 
mil: 
CDo; 
CP. 
Cle. 
CLalalpha; 
CLa; 
CLa; 
CLae; 
CMo: 
CMa: 
CMa: 
CMda: 
CMde; 
CY: 
CY¥dr 
CLr; 
Clap. 
CLda; 
CEb: 
CLdar: 
CNda: 
CNb;: 
CNp; 
CNr; 
CNadr; 
def_rud; 
def_ail; 
def_elev; 


} FLIGHTSPECS; 


/*\:jet O:prop*/ 

/*wing span */ 

/*wing area */ 

/*mean aerodynamic cord*/ 

/*mass ei 

/*x axis inertial term*/ 

/*y axis inertial term*/ 

/*Z axis inertial term*/ 

/*cross inertial term*/ 

/*used for props*/ 

/*max afterburner thrust*/ 

/*max non afterbumer thrust */ 

/*reference drag coefficient oy 

/*W-force drag coefficient due to aoa ay 
/*reference lift coefficient “i 

/*delta alpha */ 

/*lift due to pitch moment*/ 

/*lift due to angle of attack */ 

/*lift due to elevator deflection*/ 

/*pitch moment coefficient at 0 alpha = */ 
/*pitch moment coefficient due to pitch rate */ 
/*pitch moment coefficient due to aoa oh 
/*pitch moment coefficient due to delta aoa */ 
/*pitch moment coefficient from elevator motion*/ 
/*Y-force coefficent due to side slip angle */ 
/*Y-force coefficent due to change in yaw rate */ 
/*roll moment coefficient due to yaw rate */ 
/*roll moment coefficient due to roll rate */ 
/*roll moment coefficient due to change 1n aul */ 
/*roll moment coefficient due to side slip */ 
/*roll moment coefficient due to change in yaw */ 
/*aoa change effect on yaw moment + 
/*sideslip change effect on yaw +) 

/*roll rate effect on yaw moment +} 

/*yaw rate effect on yaw moment a 

/*yaw acceleration effect on yaw moment  */ 
/*radder deflection +/- in radians 7 
/*aileron deflection +/- inradians  */ 
/*elevator deflection +/- in radians oy, 





Figure C.3: Aircraft Specification Data Structure 


C. INITIALIZATION OF THE AIRCRAFT DATA STRUCTURE 


To give the system user a convenient method of assigning aircraft specification data 
to each aircraft in NPSNET, a data input file (aero.dat), together with “read_acftdat.c”’, is 


used. One simply enters the data for a particular type of aircraft in the data input file and 
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specifies, in the type aircraft field, a particular vehicle type. The vehicle types have been 
previously defined in NPSNET for distinguishing between types of vehicles. If an aircraft 
in NPSNET has a type number that does not match a type number entered in the data file 
then a default type is assigned by way of the “read_acftdat.c” program. This default type 
provides the “easiest” aircraft to pilot. An unlimited number of types can be entered 
through this data file. Users should insure that the number of these records are specified at 


the top of the file (Figure A3.4). 


flabel 
/# of aircraft types 


/record type/ 

/match to type of aircraft/ 
/jet or prop jet:1 prop:0/ 
[o/ 


/max thrust/ 
/muil thrust/ 
/CDo /CDa 
/CLo /CLa /CLq /CLda /CLde 
0.0 -3.6 -0.38 -1.1 -0.5 /CMo /CMg /CMa /CMda /CMde 
-0.98 0.17 /CYbD/CYdr 
-0.12 -0.26 0.14 0.08 -0.105 /CLb /CLp /CLr /CLda /CLdr 
0.250022 -0.35 0.06 0.032 /CNb /CNp /CNr /CNda /CNdr 
0.2618 -0.5236 0.5236 /deflection of rud, ail + stab 





Figure C.4: Aircraft Data Record 
“Read_acftdat.c” contains the code for linking the aircraft structures together. It takes 
the entire vehicle array, determines which vehicles are aircraft, matches the type records in 


the data file with the specified vehicle types, initializes the two data structures (above), and 


places an updated vehicle structure back into NPSNET. 


D. COCKPIT DISPLAY 


In order to have an appropriate interface between the user and the aircraft simulation, 
a rudimentary cockpit was devised and inserted in the code at the end of the procedure 


“graphicsproc()” found in jeep.c. 
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The projection of the cockpit is ortho2 and the following code was added: 


prefposition(0, XMAXSCREEN,0, YMAXSCREEN 
winset(jeepwin): 
acft_cockpit(vehpos); 


For control of the radar sweep the following global variables were added: 


float deltatime; /*frames per second*/ 

float radar_pos: /*keeps track of radar look angle*/ 
float rdr_posx; /*radar sweep X axis position*/ 

int reverse_dir; /*radar scan direction*/ 


The following definitions were placed 1n the “aero.h” file: 


LEFIYSCOPE 43.5 

RIGHTSCOPE 59.0 

TOPSCOPE 14.0 

BOTTOMSCOPE 3.5 

SCANRATE 72.0 

SCANWIDTH 120.0 

SCOPEWIDTH (RIGHTSCOPE-LEFTSCOPE) 
HALFSCANWIDTH (0.5*SCANWIDTH) 
SCOPERATE (SCANWIDTH/SCOPEWIDTH) 
CENTERSCOPE (LEFTSCOPE+(SCOPEWIDTH/2.0)) 
COMPASSRADIUS 5.0 

TEN (COMPASSRADIUS-0.5) 

THIRTY (COMPASSRADIUS- 1.0) 


E. INPUT DEVICE MODIFICATION FOR AIRCRAFT CONTOL 


1. Spaceball 
Control input for an aircraft is limited to movement of the Spaceball along the 
forward-aft axis and the left-right axis. Therefore, the method in which the Spaceball data 
is read was modified. To accomplish this, the procedure, newsbvals() in jeep.c, was 
separated by an ‘if’ statement into two cases: (1) an aircraft and (2) not an aircraft. From 
the Spaceball data array, values [0] and [2] are converted into aileron and elevator position 


data. The range of values is limited to -1.0 and 1.0. 


2. Keyboard 


The keyboard has been modified to add rudder, elevator trim and thrust input. The 
“>" and “<“ keys control the rudder input. The rudder can be moved left and right in .05 


radian increments, and maintains a value between -1.0 and 1.0. 
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Elevator trum is controlled by the ‘e’/’E’ and ‘r’/’R’ keys. ‘e’/’E’ trims the 
elevator surface up and ‘r’/’R’ trims the elevator surface down. Thrust ts controlled by the 
‘s'/'S’,a A’ ,and ‘d’/’D’ keys, and has a value between 0.0 to 1.0. The ‘s’/’S’ key moves 
the throttle higher in 0.05 unit increments, the ‘d’/’D’ key moves the throttle higher in 0.3 


unit increments, and the ‘a’/’A’ key moves the throttle down in 0.05 unit increments. 


F. INTEGRATION OF THE AERODYNAMIC MODEL 


The aerodynamic model used can be found in Appendix A and 1s incorporated as a 
separate compilation file “aero.c’’. Within “jeep.c” a test is made as to whether the vehicle 
being controlled is a ground or aur vehicle. If it is an air vehicle, program control moves to 
procedure movetheaircraft() found in “jeepmot.c”, where the procedure call, aero_model(), 
is made. The aerodynamic model functions exactly as described in the prototype simulator. 
After retuming from procedure aero_model(), position upd: °s are calculated and all the 
variables used by NPSNET are assigned new values. 

One aspect of NPSNET is the ability to switch from one vehicle to another through the 
planar map display. Since the vehicle starts out with an initial heading and airspeed, it is 
unpoitant to “remember” this information when assuming control of the vehicle. This is 
accomplished by initializing the aircraft state structure (Figure C.2) with airspeed and 
attitude data within the procedure switchveh() in file “dogsncats.c’’. It is important to 
remember at this point that the aerodynamic model will convert all zero airspeed values to 
10 ft/sec. If the desire is to have an immediately flyable plane at the vehicle switch, then 


the airspeed should be initialized to at least 100 ft/sec. 


G. INTEGRATION OF ORIENTATION MODEL 

The orientation model can be found in file aero.c in two separate places. The first place 
is at the end of procedure aero_model(), where it is used to determine pitch, roll and yaw 
data for the driven aircraft. The second place is in the same file but in a separate procedure, 


rotate_translate_acft(). The procedure rotate_translate_acft() can be used _ in future 
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NPSNET improvements when a quaternion rotation model is desired for updating the 


orientation and positioning of all non-driven vehicles in NPSNET (Figure C.5). 


rotate translate acft(ACFIPTR P, float dtime) /*framerate*/ 
{ 

Matrix A_MATRIX; 

float vel wor/[3]; 


/xxxkkkincrementally rotate quaternion*®KKRRRKAKKKKKKKKKEKEKE / 
increment quat (P—->ang vel [0], Po>eng vel (il) ses ang vcl[2]; 
P—>quat,dtime) ; 


/*****update rotation matrix Erom quat@pnrom 9 9" 
rotation matrix from quat (P=>quat, P=>himeeu 7, 


/****kextract euler angles From Mat rix* RR AKKKRKKKRKKRKKKEK KK / 
euler angles from_matrix(&P—>sroll,&P—->croll, 
&P=>sptch, &P--epeen, 
&P->shdg, &P->chdg, P->h_matrix) ; 
P->euler angles[1] = (fasin((P->sptch) )); 
P->euler angles[0] = (fasin((P->sroll))); 
if (P->croll < 0.0) 
if (P->sroll < 0.0) P->euler angles[(0) = (-PI - P->euler angles([0])); 
else P->euler angles([(0] = (PI - P->euler angles[0]); 7 


P->euler angles{Z] =F (fasin( (P=>shdg) 7; 
if (P->shdq < 7070) 


{ 
if (P->chdg < 0.0) P-Seuler angles([2) =P) — euler angles |, 


else P->euler angles[2] = TWOPI + P=>euloreanglos(2); 

) 

else 

if (P->chdg < 0.0) P=>euler angles(2))= PI = P->ceulecreaAngi ee, 


[xktkKRKKKCalcUbAtG New world velocity in esque coords <* <4 75 ere 
transpose _matrix(P->h_matrix,A MATRIX); 

post_mult_matrix_by vector(A_MATRIX,P->lin vel,vel wor); 
P->vel_world[0] = vel_wor[1] * FT_TO METERS; - 

P->vel world[1] -vel wor[{2] * FT TO METERS; 

P->vel world[2] = -vel_wor[0] * FT TO METERS; 


/***kkcalculate new world position 18 Sgi Coordg* *kXKKKKKKKKKKKEE / 
P->h_matrix(3] (0) += Ps>vel world[0))4edeine; 

P->h_matrix(3) [1] += P->vel_world[1] * dtime; 

P->h_matrix(3])[(2]) += P->vel world[2] * dtime; 

P=Shomatrix| 3 lS e=al Ore 


[****kkkkcalculate world lookat POLNK RAK KK KKK KKK KK KKK K / 
P->lookatpt[{0] = P->pos[0] + (cosptch * coshdg) ; 
P->lookatpt[1] = P->pos[1] + sinptch; 

P-—>lookatpt [2] P=—>pos(2]) + (eespteche ssanhda); 





Figure C.5: Procedure rotate_translate_acft() 
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